#include "GF.h"

///////////////////////////////////////////////////////////////////////////////
// V = GF_BIN_PRIMPOLY(M)
// Return a primitive polynomial of degree m over a binary extension field 2^m.
// Performs a table lookup to POLY_TABLE_PATH.
//
// Input:
//
//  m       Return a polynomial of degree m from GF(2^m).  2 <= m <- 10,000.
//
// Output:
//
//  v       Vector representation of primitive polynomial (in order of 
//          increasing powers of x).  Vector is in the form of an integer 
//          array.
///////////////////////////////////////////////////////////////////////////////
Vector * gf_bin_primpoly(int m) {

    // Check input
    if (m < 2 || m > 10000) {
        printf("Error: invalid primitive polynomial order m. Polynomial order "
               "must be between 2 and 10000.\n");
        return NULL;
    }

    // Try to open the file at path POLY_TABLE_PATH
    FILE * fp;      // File pointer
    char line[80];  // Line buffer
    int i;          
   
    //Init line
    memset(line, '\0', 80);
 
    if (( fp = fopen(POLY_TABLE_PATH, "r")) == NULL) {
        printf("Error: Could not open the polynomial table file at %s\n",
                                                           POLY_TABLE_PATH);
        return NULL;
    }

    // Read through lines until the line containing the desired polynomial
    // order is reached.
    for (i = 0; i < m-1; i++) {
        if (fgets(line, 80, fp) == NULL) {
            break;
        }
    }

    // Close the polynomial table file
    fclose(fp);

    // Check to make sure that i = m-1.  If it does not, there was an error
    // in file reading
    if (i != m-1) {
        printf("Error reading polynomial table file at %s: Could not read past "
               "line %i\n",POLY_TABLE_PATH,i);
    }

    // Count the number of commas in line.  This indicates whether the 
    // polynomial has three terms or five.
    int comma_count;
    comma_count = 0;
    for (i = 0; i < 80; i++) {
        if (line[i] == ',') comma_count++;
    }

    // Allocate a vector of sufficient size to return the polynomial vector
    Vector *poly;
    poly = allocate_vector((unsigned long)(m+1));

    // Allocate a vector to hold the polynomial exponents (either 3 or 5)
    Vector *exps;
    if (comma_count == 1) {
        exps = allocate_vector(3);
        for (i = 0; i < 3; i++) exps->data[i] = 0;
    }
    else {
        exps = allocate_vector(5);
	for (i = 0; i < 5; i++) exps->data[i] = 0;
    }

    // Populate the exps vector by separating line on commas
    char *exp_value;
    i = 0;
    exp_value = strtok(line,",");
    while (exp_value != NULL) {
        exps->data[i] = atoi(exp_value);
        i++;
        exp_value = strtok(NULL,",");
    }
    exps->data[i] = 0;

    // Populate the poly vector
    int j;  // Index into the exps array
    j = 0;
    for (i = m; i >= 0; i--) {
        if (i == exps->data[j]) {
            poly->data[i] = 1;
            j++;
        }
        else {
            poly->data[i] = 0;
        }
    } 

    // Free the exponents vector
    free_vector(exps);

    // Return the polynomial vector
    return poly;
}

///////////////////////////////////////////////////////////////////////////////
// V = GF_VECTOR(M, PRIM_POLY)  Construct the field of GF(2^m) in vector form 
// using the primitive polynomial stored in vector form in prim_poly
//
// Input:
//e
//    m          Construct the field GF(2^m)
//
//    prim_poly  Primitive m'th order polynomial in GF(2^m).  Use a root of
//               this polynomial to construct the vector space.
//
//  Output:
//
//    v          Vector space representation of GF(2^m) (in order of 
//               increasing powers of x.
//    
///////////////////////////////////////////////////////////////////////////////

Matrix * gf_vector(int m, Vector *prim_poly) {

    Vector *g, *p;
    Matrix *V;
    unsigned long i,j,N;

    // Populate g
    g = allocate_vector(prim_poly->len);
    for (i = 0; i < prim_poly->len - 1; i++) {
        g->data[i] = prim_poly->data[i];
    }
    g->data[i] = 0;

    // Initialize p
    p = allocate_vector(g->len);
    for (i = 0; i < g->len; i++) {
        p->data[i] = 0;
    }
    p->data[0] = 1;

    // Initialize V
    N = ((unsigned long)1) << m;
    V = allocate_matrix(N, g->len - 1);
    for (i = 0; i < N; i++) {
        for (j = 0; j < g->len - 1; j++) {
            V->data[i][j] = 0;
        }
    }

    // Generate the vector space
    for (i = 0; i < N-1; i++) {

        // Shift p to the right
        for (j = p->len - 1; j > 0; j--) {
            p->data[j] = p->data[j-1];
        }
        p->data[0] = 0;

        if (p->data[p->len - 1] == 1) {
            p->data[p->len - 1] = 0;
            for (j = 0; j < p->len; j++) {
                p->data[j] = p->data[j] ^ g->data[j];
            }
        }
        for (j = 0; j < p->len - 1; j++) {
            V->data[i][j] = p->data[j];
        }
    }
   
    // Free g and p
    free_vector(g);
    free_vector(p);

    return V;
}

///////////////////////////////////////////////////////////////////////////////
// A = GF_ADD_ONE(M,V) Build an add-one table for the GF with vector space 
// stored in v.
//
//  Input
//         v   The vector space of the Galois Field for which the add-one table
//             is to be constructed.
//
//  Output
//         A   Add-one table for this Galois Field.  Format:
//
//             log_a(a^q)    log_a(a^q + 1)
//
//             An entry of -1 indicates a value for which the log is not 
//             defined.
///////////////////////////////////////////////////////////////////////////////

Matrix * gf_add_one(int m, Matrix * v) {
   
    unsigned long i,j,k;
    Matrix * A;
    Vector * t;
    int match;

    //Allocate and initialize the add-one table 
    A = allocate_matrix(v->rows , (unsigned long)2);

    // Populate the table
    t = allocate_vector(v->cols);

    for (i = 0; i < v->rows - 2; i++) {
        A->data[i][0] = (int)i + 1; //FIXME this could overflow!
        
        // Copy the i'th row of v into t
        for (k = 0; k < v->cols; k++) {
            t->data[k] = v->data[i][k];
        }
        
        //add one
        t->data[0] = ((int)1) ^ (t->data[0]);

        // Match to entry in v
        for (j = 0; j < v->rows; j++) {
            match = 1;
            for (k = 0; k < v->cols && match; k++) {
                match = (t->data[k] == v->data[j][k]);
            }
            if (match) {
                A->data[i][1] = (int)j + 1; //FIXME this could overflow!
            }
        }
    }
  
    // Free the temporary vector t
    free_vector(t);
 
    // Set the zero logs
    A->data[A->rows - 2][0] = 0;
    A->data[A->rows - 2][1] = -1;
    A->data[A->rows - 1][0] = -1;
    A->data[A->rows - 1][1] = 0;
 
    return A;
}

///////////////////////////////////////////////////////////////////////////////
// Z = ELEMENT_REDUCE(X)    Reduce a sum of GF elements into a single GF
//                          element
//
//    Inputs
//
//            x    A vector of powers of GF elements to be reduced.  No
//                 particular order is required.  All elements must be
//                 nonnegative integers between 0 and 2^m - 1.
//
//    Outputs
//
//            z    A single-element vector which is a power of alpha.  The 
//                 formulation is:
//                      a^z = a^x1 + a^x2 + ... + a^xn
///////////////////////////////////////////////////////////////////////////////

int elementReduce(int m, Matrix * add_one, Vector * x) {
    int z;
    int i,j;

    // Reduce all values mod 2^m - 1
    int modval = ((1<<m) - 1); // Use bitshift to compute 2^m - 1
    for (j = 0; j< x->len; j++){
       x->data[j] = x->data[j] % modval;
    }

    // If the vector is length 1, return its data
    if (x->len == 1) {
        return x->data[0];
    }

    // Remove any zeros from the vector
    int zero_count;
    zero_count = 0;
    for (j = 0; j < x->len; j++) {
        if (x->data[j] == -1) {
            zero_count++;            
        }
        else {
            x->data[j-zero_count] = x->data[j];
        }
    }

    // If any zeros were removed, reallocate and adjust len
    int *newdata;
    if (zero_count > 0) {
        newdata = (int *)malloc((x->len - zero_count) * sizeof(int));
        memcpy(newdata, x->data, (x->len - zero_count) * sizeof(int));
        free(x->data);
        x->data = newdata;
        //free(newdata);
        x->len = x->len - zero_count;
    }

    // Sort the addends into descending order
    qsort (x->data, x->len, sizeof(int), (void *)compareFunction);

    // Combine any two terms which are the same power of alpha
    int *newdata2;
    if (x->len == 1) {
        return x->data[0];
    }
    else if (x->len > 1) {
        if (x->len == 2) {
            if (x->data[0] == x->data[1]) {
                return -1;
            }
        }
        else {
            newdata = (int *)malloc((x->len) * sizeof(int));
            i = 0;
            j = 0;
            while (i < x->len) {
                if (i < x->len - 1 && x->data[i] == x->data[i+1]) {
                    i += 2;
                }
                else {
                    newdata[j] = x->data[i];
                    j++;
                    i++;
                }
            }
            newdata2 = (int*)malloc(j * sizeof(int));
            for (i = 0; i < j; i++) {
                newdata2[i] = newdata[i];
            }
            free(x->data);
            free(newdata);
            x->data = newdata2;
            x->len = j;
        }
    }

    // If all elements have canceled, return -1
    if (x->len == 0) {
        return -1;
    }

    // Loop Invariant: x contains a^n1 + a^n2 + ...
    int diff, added;
    for(j=0; j < x->len-1; j++) {
	diff = x->data[j] - x->data[j+1];
	added = add_one->data[diff-1][1];
	x->data[j+1]=elementMult(m, add_one, added, x->data[j+1]);
    }

    if (x->len < 1)
       z = -1;
    else
       z = x->data[x->len-1];

    return z;
}

// Using a model from http://www.cs.utah.edu/dept/old/texinfo/glibc-manual-0.02/library_8.html
// I'm not entirely sure if this will work because the original qsort was setup to sort in
// ascending order... I'm trying to trick it by swapping the comparison result.
// Seemed to work when tested using test code
int compareFunction (const int *a, const int *b){
	int temp = *a - *b;

	if (temp < 0)
		return 1;
	else if (temp > 0)
		return -1;
	else
		return 0;
}

/////////////////////////////////////////////////////////////////////////////
//Z = ELEMENTMULT(X,Y) Multiply two Galois Field elements.
//
//   Inputs:
//
//      X,Y     Exponents of the two multiplicands.  The evaluated product
//               is expected to be of the form (a^X)(a^Y)
//
//   Outputs:
//
//       Z       The exponent of the product of X and Y.
////////////////////////////////////////////////////////////////////////////

int elementMult(int m, Matrix * add_one, int x, int y) {
    int z;
    int power;

    if (x == -1 || y == -1) {
        return -1;
    }
    else{
	// Calculate the modulo divisor
	power = ((1<<m) - 1);
	
	// Perform multiplication modulo divisor
	z = (x+y)%power;
    }

    // Return result
    return z;
}

///////////////////////////////////////////////////////////////////////////
// Y ISZERO(E)  Determine whether the GF element e is zero.
//
//   Inputs
//   
//       e   Element to be tested
//
//   Outputs
//
//      y    True if the element is zero (an empty array), false if not.
//
//   Notes
//       This function is simply a wrapper function for isempty().  It is
//       used solely to increase code readability.
//////////////////////////////////////////////////////////////////////////
int isZero(int m, Matrix * add_one, Vector * e){
	int intE;
	intE = elementReduce(m,add_one, e);
	return (intE < 0);
}

//////////////////////////////////////////////////////////////////////	
//ELEMENTADD(X,Y) ADD two Galois Field elements.
//
//   Inputs:
//
//       X,Y     Exponents of the two addends.  The evaluated product
//               is expected to be of the form (a^X) + (a^Y)
//
//   Outputs:
//
//       Z       The exponent of the sum of X and Y.
/////////////////////////////////////////////////////////////////////

int elementAdd(int m, Matrix * add_one, int x, int y){
	int z;
	Vector *reduceVect, *xVect, *yVect;

	// Make x and y into a vector
	reduceVect = allocate_vector(2);
	xVect = allocate_vector(1);
	yVect = allocate_vector(1);
	xVect->data[0] = x;
	yVect->data[0] = y;
	reduceVect->data[0] = x;
	reduceVect->data[1] = y;

	if (isZero(m, add_one, xVect)){
		z = y;
	}
	else if(isZero(m, add_one, yVect)){
		z = x;
	}
	else{
		// Perform element reduction to calculate z
	   z = elementReduce(m, add_one, reduceVect);
	}

	// Free Vectors
	free_vector(reduceVect);
	free_vector(xVect);
	free_vector(yVect);

	// Return result of elementAdd
	return z;	
}

////////////////////////////////////////////////////////////////////////////////
//ELEMENTDIV(N,D) Divide two Galois Field elements.
//
//   Inputs:
//
//       N,D     Exponents of the numerator (N) and denominator (D).  
//                The evaluated product is expected to be of the form (a^N)/(a^D)
//
//   Outputs:
//
//       Z       The exponent of the quotient of a^N / a^D.
///////////////////////////////////////////////////////////////////////////////


int elementDiv(int m, Matrix * add_one, int n, int d) { //Vector * n, Vector *d){

    int q; // Quotient (to be returned)
    
    // Check for divide-by-zero error
    if (d == -1) {
        printf("ERROR: DIVIDE BY ZERO IN ELEMENTDIV!\n");
        return 0;
    }
    else if (n == -1) {
        return -1;
    }
    else {
        if (n >= d) {
            q = n - d;
        }
        else {
            while (n < d) {
                n = n + (1<<m) - 1;
            }
            q = n - d;
        }
    }

    return q;

}



////////////////////////////////////////////////////////
// Z = ELEMENTPOW(X,N) Evaluate x^n
//
//   Inputs:
//
//       x       GF element
//       n       Power.  
//
//
//   Outputs:
//
//       Z       The quantity x^n
////////////////////////////////////////////////////////
int elementPow(int m, Matrix * add_one, int x, int n){
    int z;
    Vector * vectX;

    // Check if x is zero
    if (x == -1) {
       z = -1;
    }
    else if (n == 0) {
       z = 0;
    }
    else{
       vectX = allocate_vector(1);
       vectX->data[0]= x * n;
       z = elementReduce(m, add_one, vectX);
       free_vector(vectX);
    }

    return z;
}

